Data analysis for:
Why are telomeres the length that they are? Insight from a phylogenetic comparative analysis Workflow of the analyses produced by Dylan J. Padilla Perez, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.Library
library(AICcmodavg)
library(ape)
library(caper)
library(car)
library(coda)
library(extrafont)
library(geiger)
library(kableExtra)
library(MuMIn)
library(nlme)
library(pbapply)
library(phylopath)
library(phytools)
library(plotrix)
library(rphylopic)
library(scales)
library(shape)
library(xtable)R.version> _
> platform x86_64-apple-darwin17.0
> arch x86_64
> os darwin17.0
> system x86_64, darwin17.0
> status
> major 4
> minor 2.2
> year 2022
> month 10
> day 31
> svn rev 83211
> language R
> version.string R version 4.2.2 (2022-10-31)
> nickname Innocent and TrustingsessionInfo()> R version 4.2.2 (2022-10-31)
> Platform: x86_64-apple-darwin17.0 (64-bit)
> Running under: macOS Catalina 10.15.7
>
> Matrix products: default
> BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rmarkdown_2.21 xtable_1.8-4 shape_1.4.6 scales_1.2.1
> [5] rphylopic_1.2.0 plotrix_3.8-2 phylopath_1.1.3 pbapply_1.7-2
> [9] nlme_3.1-162 MuMIn_1.47.5 kableExtra_1.3.4 geiger_2.0.11
> [13] phytools_2.0-3 maps_3.4.1 extrafont_0.19 coda_0.19-4
> [17] car_3.1-2 carData_3.0-5 caper_1.0.1 mvtnorm_1.1-3
> [21] MASS_7.3-59 ape_5.7-1 AICcmodavg_2.3-2
>
> loaded via a namespace (and not attached):
> [1] subplex_1.8 doParallel_1.0.17 webshot_0.5.4
> [4] httr_1.4.5 numDeriv_2016.8-1.1 tools_4.2.2
> [7] bslib_0.4.2 utf8_1.2.3 R6_2.5.1
> [10] colorspace_2.1-0 tidyselect_1.2.0 mnormt_2.1.1
> [13] phangorn_2.11.1 grImport2_0.2-0 curl_5.0.0
> [16] compiler_4.2.2 extrafontdb_1.0 cli_3.6.1
> [19] rvest_1.0.3 expm_0.999-7 xml2_1.3.4
> [22] unmarked_1.2.5 sass_0.4.5 quadprog_1.5-8
> [25] systemfonts_1.0.4 stringr_1.5.0 digest_0.6.31
> [28] svglite_2.1.1 base64enc_0.1-3 jpeg_0.1-10
> [31] pkgconfig_2.0.3 htmltools_0.5.5 fastmap_1.1.1
> [34] highr_0.10 rlang_1.1.1 rstudioapi_0.14
> [37] VGAM_1.1-8 optimParallel_1.0-2 farver_2.1.1
> [40] jquerylib_0.1.4 generics_0.1.3 combinat_0.0-8
> [43] jsonlite_1.8.4 dplyr_1.1.2 magrittr_2.0.3
> [46] Matrix_1.5-4 Rcpp_1.0.10 munsell_0.5.0
> [49] fansi_1.0.4 abind_1.4-5 lifecycle_1.0.3
> [52] scatterplot3d_0.3-43 stringi_1.7.12 yaml_2.3.7
> [55] clusterGeneration_1.3.7 plyr_1.8.8 grid_4.2.2
> [58] parallel_4.2.2 lattice_0.21-8 splines_4.2.2
> [61] knitr_1.42 pillar_1.9.0 igraph_1.4.2
> [64] codetools_0.2-19 stats4_4.2.2 fastmatch_1.1-3
> [67] XML_3.99-0.14 glue_1.6.2 evaluate_0.20
> [70] deSolve_1.35 vctrs_0.6.2 png_0.1-8
> [73] foreach_1.5.2 Rttf2pt1_1.3.12 gtable_0.3.3
> [76] cachem_1.0.8 ggplot2_3.4.2 xfun_0.39
> [79] rsvg_2.4.0 survival_3.5-5 viridisLite_0.4.2
> [82] tibble_3.2.1 iterators_1.0.14set.seed(80)
## Dataset
data <- read.csv("new_data_2024.csv")
str(data)> 'data.frame': 159 obs. of 18 variables:
> $ Common.name : chr "Adelie penguin" "Alpine swift " "American bald eagle " "American flamengo" ...
> $ Scientific_name : chr "Pygoscelis_adeliae" "Tachymarptis_melba" "Haliaeetus_leucocephalus" "Phoenicopterus_ruber" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Aves" "Aves" "Aves" "Aves" ...
> $ Order : chr "Sphenisciformes" "Apodiformes" "Accipitriformes" "Phoenicopteriformes" ...
> $ Family : chr "Spheniscidae" "Apodidae" "Accipitridae" "Phoenicopteridae" ...
> $ Genus : chr "Pygoscelis" "Tachymarptis" "Haliaeetus" "Phoenicopterus" ...
> $ Species : chr "adeliae" "melba" "leucocephalus" "ruber" ...
> $ Endo_ectotherm : chr "endo" "endo" "endo" "endo" ...
> $ Adult_mass_grams : num 5020 103 3175 3066 117 ...
> $ log_mass : num 3.7 2.01 3.5 3.49 2.07 ...
> $ Lifespan_years : num 20 26 48 44 17 25 16 28.2 29.9 43.7 ...
> $ Loglifespan : num 1.3 1.41 1.68 1.64 1.23 ...
> $ Average_Telomere_Length_kb: num 7.35 18.78 2.1 18.31 30.2 ...
> $ logTL : num 0.866 1.274 0.322 1.263 1.48 ...
> $ Telomerase_activity : int NA NA NA NA NA NA NA NA NA NA ...head(data)> Common.name Scientific_name Domain Kingdom Phylum
> 1 Adelie penguin Pygoscelis_adeliae Eukaryota Metazoa Chordata
> 2 Alpine swift Tachymarptis_melba Eukaryota Metazoa Chordata
> 3 American bald eagle Haliaeetus_leucocephalus Eukaryota Metazoa Chordata
> 4 American flamengo Phoenicopterus_ruber Eukaryota Metazoa Chordata
> 5 American kestrel Falco_sparverius Eukaryota Metazoa Chordata
> 6 Andouin's gull Ichthyaetus_audouinii Eukaryota Metazoa Chordata
> Class Order Family Genus Species
> 1 Aves Sphenisciformes Spheniscidae Pygoscelis adeliae
> 2 Aves Apodiformes Apodidae Tachymarptis melba
> 3 Aves Accipitriformes Accipitridae Haliaeetus leucocephalus
> 4 Aves Phoenicopteriformes Phoenicopteridae Phoenicopterus ruber
> 5 Aves Falconiformes Falconidae Falco sparverius
> 6 Aves Charadriiformes Laridae Ichthyaetus audouinii
> Endo_ectotherm Adult_mass_grams log_mass Lifespan_years Loglifespan
> 1 endo 5020.0 3.700704 20 1.301030
> 2 endo 102.7 2.011570 26 1.414973
> 3 endo 3175.0 3.501744 48 1.681241
> 4 endo 3066.0 3.486572 44 1.643453
> 5 endo 117.0 2.068186 17 1.230449
> 6 endo 535.0 2.728354 25 1.397940
> Average_Telomere_Length_kb logTL Telomerase_activity
> 1 7.350 0.8662873 NA
> 2 18.783 1.2737650 NA
> 3 2.100 0.3222193 NA
> 4 18.310 1.2626883 NA
> 5 30.200 1.4800069 NA
> 6 26.550 1.4240645 NAunique(data$Class)> [1] "Aves" "Fish" "Mammalia" "Reptilia"dat <- data
names(dat)> [1] "Common.name" "Scientific_name"
> [3] "Domain" "Kingdom"
> [5] "Phylum" "Class"
> [7] "Order" "Family"
> [9] "Genus" "Species"
> [11] "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years"
> [15] "Loglifespan" "Average_Telomere_Length_kb"
> [17] "logTL" "Telomerase_activity"dat$log_mass <- log1p(dat$Adult_mass_grams)
dat$log_mass> [1] 8.5213844 4.6415021 8.0633778 8.0284552 4.7706846 6.2841342
> [7] 2.9601051 7.4312997 6.0602907 8.0811658 5.7620514 4.6463121
> [13] 6.6603190 4.7957905 3.9796817 6.1758673 7.4809922 4.3782696
> [19] 6.4393504 2.5649494 7.2449415 2.9391619 9.0669318 8.9142228
> [25] 5.5093883 7.4679423 3.8199077 6.6427475 2.7725887 5.2729996
> [31] 6.8123451 8.3696208 6.9517722 5.9162021 7.2174434 2.8903718
> [37] 6.9255952 4.9199809 2.6173958 3.0540012 8.2158178 6.8892246
> [43] 6.8721281 4.8751973 2.9957323 8.8604992 2.5877640 8.7079788
> [49] 2.5649494 6.2841342 6.3985949 2.8273136 8.3022658 0.5306283
> [55] 8.0067008 7.7836406 7.3783837 9.2687036 0.4054651 1.0296194
> [61] 1.5475625 1.0986123 0.8329091 2.5877640 10.8742854 8.7949764
> [67] 8.0712185 6.5191473 9.4880479 8.4734503 9.4121377 8.0394799
> [73] 8.0067008 8.2942996 7.2717037 8.5183925 7.4091364 10.5966597
> [79] 5.4205350 18.4206808 17.1654147 13.5923683 11.2897944 11.5253579
> [85] 11.0186455 12.2060776 12.9808021 11.7752974 9.7981826 9.9523253
> [91] 13.5278298 12.2783980 12.1007177 11.0509059 10.5966597 13.0710722
> [97] 8.2689882 8.3723986 11.6927522 9.4727816 2.3978953 2.6026897
> [103] 5.8607862 6.7345917 8.0067008 6.2166061 1.9878743 2.3025851
> [109] 8.1889669 8.3371091 7.9013774 7.4960973 4.6151205 3.7135721
> [115] 4.0775374 14.5925397 12.8584004 12.6115411 12.4292202 10.2576945
> [121] 8.7404967 11.0354701 8.8913740 6.8308742 10.5947830 11.0740483
> [127] 10.7140844 11.8482756 7.8461988 8.9763886 9.0162701 8.7584124
> [133] 15.3841267 14.9717629 9.9159595 10.9151066 6.2803958 7.0264268
> [139] 3.0680529 3.5835189 5.7071103 3.0680529 5.2922993 5.3033049
> [145] 9.5888453 11.9183972 1.7917595 3.0445224 9.2679487 5.7868974
> [151] 7.2647302 5.9210421 6.4669220 8.9683962 3.9926809 4.6634391
> [157] 7.2174434 8.2942996 4.5747110dat[dat$Average_Telomere_Length_kb > 200, ]> Common.name Scientific_name Domain Kingdom Phylum Class
> 107 Iberian shrew Sorex_granarius Eukaryota Metazoa Chordata Mammalia
> Order Family Genus Species Endo_ectotherm Adult_mass_grams
> 107 Eulipotyphla Soricidae Sorex granarius endo 6.3
> log_mass Lifespan_years Loglifespan Average_Telomere_Length_kb logTL
> 107 1.987874 3 0.4771213 213 2.32838
> Telomerase_activity
> 107 1#str(dat)
dat <- dat[!dat$Scientific_name == "Sorex_granarius", ] ## potential outlier
rownames(dat) <- dat$Scientific_name
head(dat)> Common.name Scientific_name
> Pygoscelis_adeliae Adelie penguin Pygoscelis_adeliae
> Tachymarptis_melba Alpine swift Tachymarptis_melba
> Haliaeetus_leucocephalus American bald eagle Haliaeetus_leucocephalus
> Phoenicopterus_ruber American flamengo Phoenicopterus_ruber
> Falco_sparverius American kestrel Falco_sparverius
> Ichthyaetus_audouinii Andouin's gull Ichthyaetus_audouinii
> Domain Kingdom Phylum Class Order
> Pygoscelis_adeliae Eukaryota Metazoa Chordata Aves Sphenisciformes
> Tachymarptis_melba Eukaryota Metazoa Chordata Aves Apodiformes
> Haliaeetus_leucocephalus Eukaryota Metazoa Chordata Aves Accipitriformes
> Phoenicopterus_ruber Eukaryota Metazoa Chordata Aves Phoenicopteriformes
> Falco_sparverius Eukaryota Metazoa Chordata Aves Falconiformes
> Ichthyaetus_audouinii Eukaryota Metazoa Chordata Aves Charadriiformes
> Family Genus Species
> Pygoscelis_adeliae Spheniscidae Pygoscelis adeliae
> Tachymarptis_melba Apodidae Tachymarptis melba
> Haliaeetus_leucocephalus Accipitridae Haliaeetus leucocephalus
> Phoenicopterus_ruber Phoenicopteridae Phoenicopterus ruber
> Falco_sparverius Falconidae Falco sparverius
> Ichthyaetus_audouinii Laridae Ichthyaetus audouinii
> Endo_ectotherm Adult_mass_grams log_mass
> Pygoscelis_adeliae endo 5020.0 8.521384
> Tachymarptis_melba endo 102.7 4.641502
> Haliaeetus_leucocephalus endo 3175.0 8.063378
> Phoenicopterus_ruber endo 3066.0 8.028455
> Falco_sparverius endo 117.0 4.770685
> Ichthyaetus_audouinii endo 535.0 6.284134
> Lifespan_years Loglifespan Average_Telomere_Length_kb
> Pygoscelis_adeliae 20 1.301030 7.350
> Tachymarptis_melba 26 1.414973 18.783
> Haliaeetus_leucocephalus 48 1.681241 2.100
> Phoenicopterus_ruber 44 1.643453 18.310
> Falco_sparverius 17 1.230449 30.200
> Ichthyaetus_audouinii 25 1.397940 26.550
> logTL Telomerase_activity
> Pygoscelis_adeliae 0.8662873 NA
> Tachymarptis_melba 1.2737650 NA
> Haliaeetus_leucocephalus 0.3222193 NA
> Phoenicopterus_ruber 1.2626883 NA
> Falco_sparverius 1.4800069 NA
> Ichthyaetus_audouinii 1.4240645 NA## Trees
full_data_tree <- read.tree("spp.tree.nwk")
is.ultrametric(full_data_tree)> [1] FALSEfull_data_tree <- force.ultrametric(full_data_tree)> ***************************************************************
> * Note: *
> * force.ultrametric does not include a formal method to *
> * ultrametricize a tree & should only be used to coerce *
> * a phylogeny that fails is.ultrametric due to rounding -- *
> * not as a substitute for formal rate-smoothing methods. *
> ***************************************************************is.ultrametric(full_data_tree)> [1] TRUEfull_data_tree>
> Phylogenetic tree with 158 tips and 157 internal nodes.
>
> Tip labels:
> Raja_montagui, Torpedo_marmorata, Urolophus_paucimaculatus, Squalus_acanthias, Mustelus_asterias, Callorhinchus_milii, ...
> Node labels:
> , '36', '37', '13', '14', '25', ...
>
> Rooted; includes branch lengths.## Full_tree
full_data_tree$tip.label> [1] "Raja_montagui" "Torpedo_marmorata"
> [3] "Urolophus_paucimaculatus" "Squalus_acanthias"
> [5] "Mustelus_asterias" "Callorhinchus_milii"
> [7] "Oncorhynchus_mykiss" "Lutjanus_argentimaculatus"
> [9] "Acanthopagrus_schlegelii" "Dicentrarchus_labrax"
> [11] "Macquaria_ambigua" "Platycephalus_bassensis"
> [13] "Trachurus_novaezelandiae" "Pseudocaranx_wrighti"
> [15] "Oreochromis_niloticus" "Oryzias_latipes"
> [17] "Fundulus_heteroclitus" "Xiphophorus_hellerii"
> [19] "Xiphophorus_couchianus" "Xiphophorus_maculatus"
> [21] "Nothobranchius_kadleci" "Nothobranchius_furzeri"
> [23] "Nothobranchius_rachovii" "Upeneichthys_vlamingii"
> [25] "Gadus_morhua" "Carassius_auratus"
> [27] "Cyprinus_carpio" "Danio_rerio"
> [29] "Anguilla_rostrata" "Didelphis_virginiana"
> [31] "Sorex_granarius" "Sorex_araneus"
> [33] "Atelerix_albiventris" "Crocuta_crocuta"
> [35] "Felis_catus" "Panthera_tigris"
> [37] "Canis_lupus" "Meles_meles"
> [39] "Ailurus_fulgens" "Zalophus_californianus"
> [41] "Ursus_maritimus" "Equus_caballus"
> [43] "Equus_grevyi" "Tapirus_indicus"
> [45] "Ceratotherium_simum" "Camelus_dromedarius"
> [47] "Tursiops_truncatus" "Balaena_mysticetus"
> [49] "Eschrichtius_robustus" "Hexaprotodon_liberiensis"
> [51] "Bos_taurus" "Capra_hircus"
> [53] "Ovis_aries" "Muntiacus_muntjak"
> [55] "Muntiacus_reevesi" "Rangifer_tarandus"
> [57] "Giraffa_camelopardalis" "Sus_scrofa"
> [59] "Pteropus_rodricensis" "Tadarida_brasiliensis"
> [61] "Myotis_lucifugus" "Sylvilagus_aquaticus"
> [63] "Oryctolagus_cuniculus" "Lepus_californicus"
> [65] "Ochotona_princeps" "Aplodontia_rufa"
> [67] "Sciurus_carolinensis" "Peromyscus_eremicus"
> [69] "Cricetulus_griseus" "Rattus_norvegicus"
> [71] "Mus_musculus" "Mus_spretus"
> [73] "Castor_canadensis" "Heterocephalus_glaber"
> [75] "Hydrochoerus_hydrochaeris" "Lemur_catta"
> [77] "Saimiri_sciureus" "Ateles_geoffroyi"
> [79] "Gorilla_gorilla" "Homo_sapiens"
> [81] "Pan_paniscus" "Pan_troglodytes"
> [83] "Pongo_pygmaeus" "Macaca_mulatta"
> [85] "Macaca_fascicularis" "Macaca_nemestrina"
> [87] "Tupaia_tana" "Tupaia_belangeri"
> [89] "Procavia_capensis" "Loxodonta_africana"
> [91] "Elephas_maximus" "Setifer_setosus"
> [93] "Galegeeska_rufescens" "Macroscelides_proboscideus"
> [95] "Chaetophractus_vellerosus" "Myrmecophaga_tridactyla"
> [97] "Choloepus_hoffmanni" "Liasis_fuscus"
> [99] "Thamnophis_sirtalis" "Lacerta_agilis"
> [101] "Zootoca_vivipara" "Alligator_mississippiensis"
> [103] "Alligator_sinensis" "Tachymarptis_melba"
> [105] "Sula_sula" "Fregata_minor"
> [107] "Pygoscelis_adeliae" "Macronectes_giganteus"
> [109] "Macronectes_halli" "Fulmarus_glacialis"
> [111] "Calonectris_diomedea" "Oceanodroma_leucorhoa"
> [113] "Diomedea_exulans" "Thalassarche_melanophris"
> [115] "Antigone_antigone" "Sterna_hirundo"
> [117] "Ichthyaetus_audouinii" "Larus_fuscus"
> [119] "Rissa_tridactyla" "Uria_lomvia"
> [121] "Cepphus_grylle" "Calidris_alpina"
> [123] "Calidris_pugnax" "Haematopus_ostralegus"
> [125] "Phoenicopterus_ruber" "Buteo_swainsoni"
> [127] "Buteo_jamaicensis" "Haliaeetus_leucocephalus"
> [129] "Accipiter_gentilis" "Gyps_fulvus"
> [131] "Amazona_amazonica" "Strigops_habroptila"
> [133] "Acrocephalus_sechellensis" "Hirundo_rustica"
> [135] "Riparia_riparia" "Tachycineta_albilinea"
> [137] "Tachycineta_bicolor" "Parus_major"
> [139] "Turdus_merula" "Taeniopygia_guttata"
> [141] "Lonchura_striata" "Chloebia_gouldiae"
> [143] "Coloeus_monedula" "Aphelocoma_ultramarina"
> [145] "Aphelocoma_coerulescens" "Passerculus_sandwichensis"
> [147] "Vireo_olivaceus" "Falco_sparverius"
> [149] "Branta_leucopsis" "Chrysolophus_pictus"
> [151] "Meleagris_gallopavo" "Gallus_gallus"
> [153] "Colinus_virginianus" "Trachemys_scripta"
> [155] "Chrysemys_picta" "Emys_orbicularis"
> [157] "Chelonia_mydas" "Neoceratodus_forsteri"check <- name.check(full_data_tree, dat)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pruned_data_tree <- drop.tip(full_data_tree, rm_phy)
pruned_dat <- subset(dat, subset = dat$Scientific_name %in% pruned_data_tree$tip, select = names(dat))
str(pruned_dat)> 'data.frame': 142 obs. of 18 variables:
> $ Common.name : chr "Adelie penguin" "Alpine swift " "American bald eagle " "American flamengo" ...
> $ Scientific_name : chr "Pygoscelis_adeliae" "Tachymarptis_melba" "Haliaeetus_leucocephalus" "Phoenicopterus_ruber" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Aves" "Aves" "Aves" "Aves" ...
> $ Order : chr "Sphenisciformes" "Apodiformes" "Accipitriformes" "Phoenicopteriformes" ...
> $ Family : chr "Spheniscidae" "Apodidae" "Accipitridae" "Phoenicopteridae" ...
> $ Genus : chr "Pygoscelis" "Tachymarptis" "Haliaeetus" "Phoenicopterus" ...
> $ Species : chr "adeliae" "melba" "leucocephalus" "ruber" ...
> $ Endo_ectotherm : chr "endo" "endo" "endo" "endo" ...
> $ Adult_mass_grams : num 5020 103 3175 3066 117 ...
> $ log_mass : num 8.52 4.64 8.06 8.03 4.77 ...
> $ Lifespan_years : num 20 26 48 44 17 25 16 28.2 29.9 28.5 ...
> $ Loglifespan : num 1.3 1.41 1.68 1.64 1.23 ...
> $ Average_Telomere_Length_kb: num 7.35 18.78 2.1 18.31 30.2 ...
> $ logTL : num 0.866 1.274 0.322 1.263 1.48 ...
> $ Telomerase_activity : int NA NA NA NA NA NA NA NA NA NA ...head(pruned_dat)> Common.name Scientific_name
> Pygoscelis_adeliae Adelie penguin Pygoscelis_adeliae
> Tachymarptis_melba Alpine swift Tachymarptis_melba
> Haliaeetus_leucocephalus American bald eagle Haliaeetus_leucocephalus
> Phoenicopterus_ruber American flamengo Phoenicopterus_ruber
> Falco_sparverius American kestrel Falco_sparverius
> Ichthyaetus_audouinii Andouin's gull Ichthyaetus_audouinii
> Domain Kingdom Phylum Class Order
> Pygoscelis_adeliae Eukaryota Metazoa Chordata Aves Sphenisciformes
> Tachymarptis_melba Eukaryota Metazoa Chordata Aves Apodiformes
> Haliaeetus_leucocephalus Eukaryota Metazoa Chordata Aves Accipitriformes
> Phoenicopterus_ruber Eukaryota Metazoa Chordata Aves Phoenicopteriformes
> Falco_sparverius Eukaryota Metazoa Chordata Aves Falconiformes
> Ichthyaetus_audouinii Eukaryota Metazoa Chordata Aves Charadriiformes
> Family Genus Species
> Pygoscelis_adeliae Spheniscidae Pygoscelis adeliae
> Tachymarptis_melba Apodidae Tachymarptis melba
> Haliaeetus_leucocephalus Accipitridae Haliaeetus leucocephalus
> Phoenicopterus_ruber Phoenicopteridae Phoenicopterus ruber
> Falco_sparverius Falconidae Falco sparverius
> Ichthyaetus_audouinii Laridae Ichthyaetus audouinii
> Endo_ectotherm Adult_mass_grams log_mass
> Pygoscelis_adeliae endo 5020.0 8.521384
> Tachymarptis_melba endo 102.7 4.641502
> Haliaeetus_leucocephalus endo 3175.0 8.063378
> Phoenicopterus_ruber endo 3066.0 8.028455
> Falco_sparverius endo 117.0 4.770685
> Ichthyaetus_audouinii endo 535.0 6.284134
> Lifespan_years Loglifespan Average_Telomere_Length_kb
> Pygoscelis_adeliae 20 1.301030 7.350
> Tachymarptis_melba 26 1.414973 18.783
> Haliaeetus_leucocephalus 48 1.681241 2.100
> Phoenicopterus_ruber 44 1.643453 18.310
> Falco_sparverius 17 1.230449 30.200
> Ichthyaetus_audouinii 25 1.397940 26.550
> logTL Telomerase_activity
> Pygoscelis_adeliae 0.8662873 NA
> Tachymarptis_melba 1.2737650 NA
> Haliaeetus_leucocephalus 0.3222193 NA
> Phoenicopterus_ruber 1.2626883 NA
> Falco_sparverius 1.4800069 NA
> Ichthyaetus_audouinii 1.4240645 NApruned_data_tree>
> Phylogenetic tree with 142 tips and 141 internal nodes.
>
> Tip labels:
> Torpedo_marmorata, Urolophus_paucimaculatus, Squalus_acanthias, Mustelus_asterias, Callorhinchus_milii, Oncorhynchus_mykiss, ...
> Node labels:
> , '36', '37', '14', '25', '286', ...
>
> Rooted; includes branch lengths.name.check(pruned_data_tree, pruned_dat)> [1] "OK"hist(pruned_dat$Lifespan_years, main = "raw variable")
hist(log1p(pruned_dat$Lifespan_years), main = "log-transformed")
names(pruned_dat)> [1] "Common.name" "Scientific_name"
> [3] "Domain" "Kingdom"
> [5] "Phylum" "Class"
> [7] "Order" "Family"
> [9] "Genus" "Species"
> [11] "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years"
> [15] "Loglifespan" "Average_Telomere_Length_kb"
> [17] "logTL" "Telomerase_activity"pruned_dat$log.lifespan <- log1p(pruned_dat$Lifespan_years)
pruned_dat$log.lifespan> [1] 3.0445224 3.2958369 3.8918203 3.8066625 2.8903718 3.2580965 2.8332133
> [8] 3.3741687 3.4307562 3.3843903 3.1267605 3.4339872 3.5263605 3.3945084
> [15] 3.7909847 2.7725887 2.6672282 3.7841896 2.7972813 3.7471484 4.1108739
> [22] 3.6109179 3.5807373 2.4510051 2.0014800 3.9512437 3.5945688 3.1354942
> [29] 3.4339872 3.4563167 2.4159138 3.1780538 2.4069451 2.9444390 3.8712010
> [36] 3.2542430 3.4011974 3.2580965 2.5726122 3.9318256 2.3978953 2.6390573
> [43] 2.5649494 3.2503745 2.8903718 3.9318256 1.7917595 2.6390573 2.3025851
> [50] 3.7376696 3.8712010 1.8718022 1.7917595 1.7917595 0.2070142 0.9162907
> [57] 1.6094379 3.2580965 2.7725887 2.0794415 2.3978953 3.0445224 2.9444390
> [64] 2.4849066 3.1780538 2.9601051 1.9459101 2.5649494 4.3307333 3.0445224
> [71] 4.4543473 2.7146947 5.3565863 4.3567088 3.7013020 3.1696856 3.1267605
> [78] 3.0819100 3.9627161 3.3809947 3.3322045 3.1863526 2.9856819 3.0445224
> [85] 3.6027768 3.7400477 3.3322045 3.8022081 3.4339872 2.9957323 3.3068867
> [92] 2.9755296 3.5553481 2.5649494 3.3672958 2.8449094 2.0281482 2.5176965
> [99] 1.4350845 2.7600099 2.5494452 2.1972246 2.6390573 2.0794415 2.2721259
> [106] 2.1860513 3.8286414 3.4657359 4.0604430 3.6243409 3.4657359 3.7376696
> [113] 4.8162412 3.8732822 3.4404181 4.0253517 4.0943446 4.2341065 4.1125119
> [120] 3.6454499 3.6532523 3.7135721 3.6888795 4.1896547 4.3894986 3.1945831
> [127] 2.7788193 3.2027464 1.9459101 3.4657359 1.5686159 1.6094379 2.5336968
> [134] 2.4932055 4.1896547 4.3567088 2.4849066 2.3978953 3.3250360 4.7957905
> [141] 3.4339872 4.1271344
Figure 1
#fonts()
#png("figure1.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure1.pdf")
par(family = "Monaco")
plot(NA, type = "n", xaxt = "n", ylab = "", xlab = "", xlim = c(0, 5), ylim = c(0, 5), axes = FALSE)
box()
text(x = 2.5, y = 5, "Body size", cex = 1.4)
text(x = 0.5, y = 0.5, "Lifespan", cex = 1.4)
text(x = 4.5, y = 0.5, "Endothermy", cex = 1.4)
text(x = 2.5, y = 2.5, "Telomere\nlength", cex = 1.4)
Arrows(2.1, 4.8, 0.7, 0.7, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")
Arrows(3, 4.8, 4.5, 0.7, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")
Arrows(1.2, 0.5, 3.5, 0.5, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")
Arrows(1, 0.7, 2.3, 2, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)
Arrows(3.9, 0.7, 2.7, 2, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)
Arrows(2.5, 4.8, 2.5, 2.8, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)
##dev.off()
Figure 2
#png("figure2.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure2.pdf")
plotTree(pruned_data_tree, fsize = 0.25, ftype = "i", mar = c(1, 1, 1, 3))
cladebox <- function(tree, node, spp, ysize, color,...){
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
h <- max(nodeHeights(tree))
parent <- tree$edge[which(tree$edge[ , 2] == node), 1]
x0 <- max(c(obj$xx[node] + obj$xx[parent])/2, obj$xx[node]-0.01*h)
x1 <- obj$x.lim[2]
dd <- getDescendants(tree, node)
y0 <- min(range(obj$yy[dd]))-0.6
y1 <- max(range(obj$yy[dd]))+0.6
cols <- colorRampPalette(c("white", color))(100)
x0 <- seq(x0, x1, length.out = 101)[1:100]
x1 <- seq(x0[1], x1, length.out = 101)[2:101]
for(j in 1:100){
polygon(c(x0[j], x1[j], x1[j], x0[j]), c(y0, y0, y1, y1), col = make.transparent(cols[j], 0.6),
border = 0)
}
spp <- gsub("_", " ", as.character(spp))
x2<- runif(100, obj$x.lim[2] + 10, obj$x.lim[2]+50)
y3 <- seq(y0, y1, 5)
for(i in spp){
x2
uuid <- get_uuid(name = i, n = 1)
img <- get_phylopic(uuid = uuid)
i <- gsub(" ", "_", i) ## This is actually not required as I think rphylopic can get spp names with "_"
nodes <- sapply(i, grep, x = tree$tip.label)
for(j in nodes){
add_phylopic_base(img = img, x = sample(x2, 1), y = j, ysize = ysize, color = color, fill = color)
}
}
}
#nodelabels(cex = 0.5)
chel.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][8]
bird.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Aves"][c(4, 10, 11, 12, 27, 28, 5)]
alg.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][2]
squ.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][4]
mammal.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Mammalia"][c(4, 5, 2, 16, 30, 37, 47, 53, 55, 49)]
ost.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Fish"][c(7, 4, 8)]
chon.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Fish"][c(19, 24)]
cladebox(pruned_data_tree, 282, spp = chel.spp, 3, color = "#e93e3a")
cladebox(pruned_data_tree, 238, spp = bird.spp, 3, color = "#759c8a")
cladebox(pruned_data_tree, 237, spp = alg.spp, 2, color = "brown")
cladebox(pruned_data_tree, 233, spp = squ.spp, 1, color = "#e7dddf")
cladebox(pruned_data_tree, 170, spp = mammal.spp, 3, color = "#dd9f40")
cladebox(pruned_data_tree, 149, spp = ost.spp, 2, color = "#a39db8")
cladebox(pruned_data_tree, 144, spp = chon.spp, 2, color = "#efe0a8")
#dev.off()
## PGLS models
names(pruned_dat)> [1] "Common.name" "Scientific_name"
> [3] "Domain" "Kingdom"
> [5] "Phylum" "Class"
> [7] "Order" "Family"
> [9] "Genus" "Species"
> [11] "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years"
> [15] "Loglifespan" "Average_Telomere_Length_kb"
> [17] "logTL" "Telomerase_activity"
> [19] "log.lifespan"pruned_dat$new.log.TL <- log(pruned_dat$Average_Telomere_Length_kb)
pruned_dat$new.log.TL> [1] 1.9947003 2.9329522 0.7419373 2.9074474 3.4078419 3.2790297 2.2586332
> [8] 2.4638532 2.6803364 2.2828925 2.4024304 2.9806186 2.5494452 2.6253930
> [15] 2.2813615 1.7101878 1.1631508 1.9586853 3.8971124 2.8534772 2.7507904
> [22] 2.9460167 2.5257286 1.8373700 3.1527360 2.6567569 2.1488510 0.2623643
> [29] 2.4680995 2.5877640 2.0918641 2.6946272 3.0568274 2.4535880 2.0893919
> [36] 3.7471484 2.0412203 2.6100698 2.6100698 2.2925348 2.7948393 2.6318888
> [43] 2.6173958 2.1587147 1.3737156 2.5257286 2.0149030 1.4011830 1.9459101
> [50] 1.5260563 1.5303947 2.7725887 1.4469190 1.4469190 1.8885837 2.0502702
> [57] 1.7917595 2.4313300 1.2354715 1.7155981 1.7011051 1.8148247 1.6601310
> [64] 2.9957323 2.5006159 1.0986123 3.2023399 3.1941732 2.5257286 1.0986123
> [71] 3.1871787 3.9120230 2.1972246 2.5649494 2.3025851 2.8903718 2.9856819
> [78] 2.3418058 2.8332133 2.5649494 2.7080502 2.6390573 2.6390573 2.8903718
> [85] 1.6094379 2.0794415 2.7080502 2.4849066 3.2695689 3.2188758 3.9120230
> [92] 2.2874715 3.4011974 3.2580965 2.1972246 2.5649494 3.5553481 3.6375862
> [99] 2.5329028 2.7080502 3.9120230 3.9120230 4.3820266 2.9957323 3.6109179
> [106] 3.6109179 2.3025851 2.4849066 2.9285235 2.4849066 2.4849066 2.8332133
> [113] 2.1972246 1.9459101 2.1972246 2.3025851 2.3025851 2.4510051 2.7343675
> [120] 2.9444390 2.7013612 2.7725887 2.7146947 2.6390573 2.7080502 1.9459101
> [127] 2.3025851 3.9120230 2.1972246 2.7725887 3.6888795 3.6888795 2.7080502
> [134] 3.1354942 2.7825391 3.4307562 2.9606231 2.9231616 3.2976870 2.9957323
> [141] 3.9120230 4.0943446model1 <- gls(new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm, correlation = corBrownian(phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(model1)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm
> Data: pruned_dat
> AIC BIC logLik
> 343.6862 370.2886 -162.8431
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.233397 1.3121680 2.4641640 0.0150
> log_mass -0.034819 0.2118296 -0.1643744 0.8697
> log.lifespan -0.138426 0.3991831 -0.3467736 0.7293
> Endo_ectothermendo -0.175596 1.7674781 -0.0993482 0.9210
> log_mass:log.lifespan 0.006902 0.0705437 0.0978350 0.9222
> log_mass:Endo_ectothermendo 0.011989 0.2326388 0.0515342 0.9590
> log.lifespan:Endo_ectothermendo 0.013390 0.4800741 0.0278912 0.9778
> log_mass:log.lifespan:Endo_ectothermendo -0.012800 0.0752810 -0.1700348 0.8652
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass -0.469
> log.lifespan -0.615 0.535
> Endo_ectothermendo -0.495 0.279 0.453
> log_mass:log.lifespan 0.456 -0.925 -0.730 -0.279
> log_mass:Endo_ectothermendo 0.427 -0.911 -0.487 -0.418 0.842
> log.lifespan:Endo_ectothermendo 0.512 -0.445 -0.832 -0.623 0.607
> log_mass:log.lifespan:Endo_ectothermendo -0.427 0.867 0.684 0.405 -0.937
> lg_:E_ lg.:E_
> log_mass
> log.lifespan
> Endo_ectothermendo
> log_mass:log.lifespan
> log_mass:Endo_ectothermendo
> log.lifespan:Endo_ectothermendo 0.562
> log_mass:log.lifespan:Endo_ectothermendo -0.922 -0.732
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.00833399 -0.09146447 0.13823694 0.30407049 0.92539247
>
> Residual standard error: 2.098608
> Degrees of freedom: 142 total; 134 residualmodel2 <- update(model1, ~ .-log_mass:log.lifespan:Endo_ectotherm)
summary(model2)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm
> Data: pruned_dat
> AIC BIC logLik
> 341.7168 365.3634 -162.8584
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.1380986 1.1821750 2.6545128 0.0089
> log_mass -0.0035886 0.1051473 -0.0341291 0.9728
> log.lifespan -0.0920198 0.2902581 -0.3170273 0.7517
> Endo_ectothermendo -0.0538850 1.6102252 -0.0334643 0.9734
> log_mass:log.lifespan -0.0043401 0.0245142 -0.1770442 0.8597
> log_mass:Endo_ectothermendo -0.0245012 0.0894885 -0.2737915 0.7847
> log.lifespan:Endo_ectothermendo -0.0463691 0.3258585 -0.1422984 0.8871
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:. lg_:E_
> log_mass -0.219
> log.lifespan -0.489 -0.160
> Endo_ectothermendo -0.389 -0.159 0.264
> log_mass:log.lifespan 0.175 -0.648 -0.349 0.317
> log_mass:Endo_ectothermendo 0.095 -0.575 0.510 -0.125 -0.165
> log.lifespan:Endo_ectothermendo 0.324 0.558 -0.666 -0.525 -0.332 -0.430
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.99962915 -0.09501665 0.14604616 0.32250953 0.93359053
>
> Residual standard error: 2.098834
> Degrees of freedom: 142 total; 135 residualmodel3 <- update(model2, ~ .-log.lifespan:Endo_ectotherm)
summary(model3)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm
> Data: pruned_dat
> AIC BIC logLik
> 339.7381 360.4289 -162.8691
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.192564 1.1144614 2.8646697 0.0048
> log_mass 0.004768 0.0869066 0.0548589 0.9563
> log.lifespan -0.119529 0.2157295 -0.5540687 0.5804
> Endo_ectothermendo -0.174090 1.3659091 -0.1274533 0.8988
> log_mass:log.lifespan -0.005499 0.0230394 -0.2386600 0.8117
> log_mass:Endo_ectothermendo -0.029972 0.0805158 -0.3722555 0.7103
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass -0.509
> log.lifespan -0.387 0.343
> Endo_ectothermendo -0.272 0.190 -0.135
> log_mass:log.lifespan 0.317 -0.591 -0.810 0.178
> log_mass:Endo_ectothermendo 0.273 -0.448 0.333 -0.456 -0.361
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.99393208 -0.09318428 0.14374430 0.32634461 0.94174318
>
> Residual standard error: 2.098991
> Degrees of freedom: 142 total; 136 residualmodel4 <- update(model3, ~ .-log_mass:log.lifespan)
summary(model4)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:Endo_ectotherm
> Data: pruned_dat
> AIC BIC logLik
> 337.7976 355.5325 -162.8988
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.276795 1.0534568 3.1105164 0.0023
> log_mass -0.007483 0.0698850 -0.1070806 0.9149
> log.lifespan -0.161249 0.1259793 -1.2799652 0.2027
> Endo_ectothermendo -0.116185 1.3395539 -0.0867343 0.9310
> log_mass:Endo_ectothermendo -0.036908 0.0748303 -0.4932183 0.6226
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct
> log_mass -0.421
> log.lifespan -0.235 -0.286
> Endo_ectothermendo -0.352 0.371 0.015
> log_mass:Endo_ectothermendo 0.438 -0.878 0.074 -0.427
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.99267537 -0.09429773 0.14396653 0.32720827 0.94298031
>
> Residual standard error: 2.099431
> Degrees of freedom: 142 total; 137 residualmodel5 <- update(model4, ~ .-log_mass:Endo_ectotherm)
summary(model5)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm
> Data: pruned_dat
> AIC BIC logLik
> 336.0495 350.8286 -163.0247
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.504551 0.9442549 3.711446 0.0003
> log_mass -0.037761 0.0333064 -1.133747 0.2589
> log.lifespan -0.156676 0.1252926 -1.250477 0.2132
> Endo_ectothermendo -0.398370 1.2079020 -0.329803 0.7420
>
> Correlation:
> (Intr) lg_mss lg.lfs
> log_mass -0.084
> log.lifespan -0.298 -0.465
> Endo_ectothermendo -0.203 -0.009 0.052
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.99465171 -0.09075665 0.13924989 0.31601039 0.93865079
>
> Residual standard error: 2.101294
> Degrees of freedom: 142 total; 138 residualmodel6 <- update(model5, ~ .-Endo_ectotherm)
summary(model6)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan
> Data: pruned_dat
> AIC BIC logLik
> 334.1614 345.9847 -163.0807
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.441364 0.9216444 3.733939 0.0003
> log_mass -0.037861 0.0331981 -1.140468 0.2561
> log.lifespan -0.154544 0.1247240 -1.239089 0.2174
>
> Correlation:
> (Intr) lg_mss
> log_mass -0.088
> log.lifespan -0.294 -0.466
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.1565570775 -0.2145557854 0.0002704793 0.1739263041 0.7765132875
>
> Residual standard error: 2.102122
> Degrees of freedom: 142 total; 139 residualmodel7 <- update(model6, ~ .-log_mass)
summary(model7)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log.lifespan
> Data: pruned_dat
> AIC BIC logLik
> 333.4839 342.3514 -163.742
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.348911 0.9190577 3.643853 0.0004
> log.lifespan -0.220773 0.1104984 -1.997977 0.0477
>
> Correlation:
> (Intr)
> log.lifespan -0.38
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.13370657 -0.20204102 -0.01077654 0.15669655 0.78439669
>
> Residual standard error: 2.111934
> Degrees of freedom: 142 total; 140 residual## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(model7), resid(model7), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals")
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(model7), col = "darkgrey")
qqline(resid(model7), col = "dodgerblue", lwd = 2)
Figure 5
#png("figure3.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure3.pdf")
par(mar = c(5, 4, 1, 4))
plot(new.log.TL ~ log.lifespan , data = pruned_dat, pch = 21, bg = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], col = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], las = 1, xlab = "Log lifespan (yrs)", ylab = "Log telomere length (kb)", type = "n")
grid(nx = NULL, ny = NULL, col = "lightgray", lwd = 1)
par(new = TRUE)
plot(new.log.TL ~ log.lifespan , data = pruned_dat, pch = 21, bg = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], col = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], las = 1, xlab = "Log lifespan (yrs)", ylab = "Log telomere length (kb)")
par(xpd = TRUE)
uuid.aves <- get_uuid(name = "Accipiter_gentilis", n = 1)
img.aves <- get_phylopic(uuid = uuid.aves)
add_phylopic_base(img = img.aves, x = 6, y = 2.4, ysize = 0.1, col = alpha("purple", 0.5), fill = alpha("purple", 0.5))
uuid.fish <- get_uuid(name = "Oreochromis niloticus", n = 1)
img.fish <- get_phylopic(uuid = uuid.fish)
add_phylopic_base(img = img.fish, x = 6, y = 2.2, ysize = 0.09, col = alpha("orange", 0.5), fill = alpha("orange", 0.5))
uuid.mammals <- get_uuid(name = "Elephas_maximus", n = 1)
img.mammals <- get_phylopic(uuid = uuid.mammals)
add_phylopic_base(img = img.mammals, x = 6, y = 2, ysize = 0.2, col = alpha("pink", 0.7), fill = alpha("pink", 0.7), horizontal = TRUE)
uuid.reptilia <- get_uuid(name = "Alligator_mississippiensis", n = 1)
img.reptilia <- get_phylopic(uuid = uuid.reptilia)
add_phylopic_base(img = img.reptilia, x = 6, y = 1.8, ysize = 0.1, col = alpha("gold", 0.5), fill = alpha("gold", 0.5))
SSX <- sum(round((pruned_dat$new.log.TL - mean(pruned_dat$new.log.TL))^2), 2)
s2 <- var(pruned_dat$new.log.TL)
n <- length(pruned_dat$new.log.TL)
x <- seq(min(pruned_dat$log.lifespan), max(pruned_dat$log.lifespan), 0.06)
m.x <- mean(round(pruned_dat$new.log.TL, 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.s
lower.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.i
##par(new = TRUE)
lines(x = x, y = (coef(model7)[1] + (coef(model7)[2]*x)), lwd = 2, col = alpha("black", 0.5))
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
Figure 4
## Phylogenetic signal (Blomberg et al.'s K statitics
telo.length <- setNames(pruned_dat$new.log.TL, rownames(pruned_dat))
telo.length> Pygoscelis_adeliae Tachymarptis_melba
> 1.9947003 2.9329522
> Haliaeetus_leucocephalus Phoenicopterus_ruber
> 0.7419373 2.9074474
> Falco_sparverius Ichthyaetus_audouinii
> 3.4078419 3.2790297
> Hirundo_rustica Branta_leucopsis
> 2.2586332 2.4638532
> Cepphus_grylle Rissa_tridactyla
> 2.6803364 2.2828925
> Turdus_merula Gallus_gallus
> 2.4024304 2.9806186
> Sterna_hirundo Calidris_alpina
> 2.5494452 2.6253930
> Haematopus_ostralegus Aphelocoma_coerulescens
> 2.2813615 1.7101878
> Chrysolophus_pictus Fregata_minor
> 1.1631508 1.9586853
> Parus_major Gyps_fulvus
> 3.8971124 2.8534772
> Strigops_habroptila Oceanodroma_leucorhoa
> 2.7507904 2.9460167
> Larus_fuscus Tachycineta_albilinea
> 2.5257286 1.8373700
> Colinus_virginianus Fulmarus_glacialis
> 3.1527360 2.6567569
> Macronectes_halli Accipiter_gentilis
> 2.1488510 0.2623643
> Amazona_amazonica Buteo_jamaicensis
> 2.4680995 2.5877640
> Vireo_olivaceus Sula_sula
> 2.0918641 2.6946272
> Riparia_riparia Passerculus_sandwichensis
> 3.0568274 2.4535880
> Macronectes_giganteus Buteo_swainsoni
> 2.0893919 3.7471484
> Uria_lomvia Aphelocoma_ultramarina
> 2.0412203 2.6100698
> Tachycineta_bicolor Diomedea_exulans
> 2.6100698 2.2925348
> Lonchura_striata Meleagris_gallopavo
> 2.7948393 2.6318888
> Taeniopygia_guttata Calonectris_diomedea
> 2.6173958 2.1587147
> Acrocephalus_sechellensis Anguilla_rostrata
> 1.3737156 2.5257286
> Oryzias_latipes Pseudocaranx_wrighti
> 2.0149030 1.4011830
> Oreochromis_niloticus Carassius_auratus
> 1.9459101 1.5260563
> Cyprinus_carpio Danio_rerio
> 1.5303947 2.7725887
> Xiphophorus_maculatus Xiphophorus_hellerii
> 1.4469190 1.4469190
> Nothobranchius_furzeri Nothobranchius_rachovii
> 1.8885837 2.0502702
> Fundulus_heteroclitus Gadus_morhua
> 1.7917595 2.4313300
> Dicentrarchus_labrax Acanthopagrus_schlegelii
> 1.2354715 1.7155981
> Upeneichthys_vlamingii Macquaria_ambigua
> 1.7011051 1.8148247
> Lutjanus_argentimaculatus Oncorhynchus_mykiss
> 1.6601310 2.9957323
> Platycephalus_bassensis Mustelus_asterias
> 2.5006159 1.0986123
> Callorhinchus_milii Urolophus_paucimaculatus
> 3.2023399 3.1941732
> Squalus_acanthias Torpedo_marmorata
> 2.5257286 1.0986123
> Neoceratodus_forsteri Setifer_setosus
> 3.1871787 3.9120230
> Balaena_mysticetus Eschrichtius_robustus
> 2.1972246 2.5649494
> Giraffa_camelopardalis Ovis_aries
> 2.3025851 2.8903718
> Rangifer_tarandus Capra_hircus
> 2.9856819 2.3418058
> Tursiops_truncatus Camelus_dromedarius
> 2.8332133 2.5649494
> Sus_scrofa Muntiacus_reevesi
> 2.7080502 2.6390573
> Muntiacus_muntjak Bos_taurus
> 2.6390573 2.8903718
> Zalophus_californianus Crocuta_crocuta
> 1.6094379 2.0794415
> Canis_lupus Ursus_maritimus
> 2.7080502 2.4849066
> Felis_catus Ailurus_fulgens
> 3.2695689 3.2188758
> Panthera_tigris Meles_meles
> 3.9120230 2.2874715
> Myotis_lucifugus Tadarida_brasiliensis
> 3.4011974 3.2580965
> Pteropus_rodricensis Chaetophractus_vellerosus
> 2.1972246 2.5649494
> Didelphis_virginiana Atelerix_albiventris
> 3.5553481 3.6375862
> Sorex_araneus Procavia_capensis
> 2.5329028 2.7080502
> Lepus_californicus Sylvilagus_aquaticus
> 3.9120230 3.9120230
> Oryctolagus_cuniculus Ochotona_princeps
> 4.3820266 2.9957323
> Macroscelides_proboscideus Galegeeska_rufescens
> 3.6109179 3.6109179
> Ceratotherium_simum Equus_grevyi
> 2.3025851 2.4849066
> Equus_caballus Tapirus_indicus
> 2.9285235 2.4849066
> Myrmecophaga_tridactyla Choloepus_hoffmanni
> 2.4849066 2.8332133
> Homo_sapiens Ateles_geoffroyi
> 2.1972246 1.9459101
> Saimiri_sciureus Pan_paniscus
> 2.1972246 2.3025851
> Pongo_pygmaeus Pan_troglodytes
> 2.3025851 2.4510051
> Gorilla_gorilla Lemur_catta
> 2.7343675 2.9444390
> Macaca_nemestrina Macaca_mulatta
> 2.7013612 2.7725887
> Macaca_fascicularis Loxodonta_africana
> 2.7146947 2.6390573
> Elephas_maximus Castor_canadensis
> 2.7080502 1.9459101
> Hydrochoerus_hydrochaeris Sciurus_carolinensis
> 2.3025851 3.9120230
> Aplodontia_rufa Heterocephalus_glaber
> 2.1972246 2.7725887
> Rattus_norvegicus Mus_musculus
> 3.6888795 3.6888795
> Tupaia_tana Tupaia_belangeri
> 2.7080502 3.1354942
> Alligator_sinensis Alligator_mississippiensis
> 2.7825391 3.4307562
> Zootoca_vivipara Lacerta_agilis
> 2.9606231 2.9231616
> Liasis_fuscus Emys_orbicularis
> 3.2976870 2.9957323
> Trachemys_scripta Chrysemys_picta
> 3.9120230 4.0943446k_tl <- phylosig(pruned_data_tree, telo.length, test = TRUE, nsim = 10000)
attributes(k_tl)> $names
> [1] "K" "P" "sim.K"
>
> $class
> [1] "phylosig"
>
> $method
> [1] "K"
>
> $test
> [1] TRUE
>
> $se
> [1] FALSEhead(k_tl$sim.K)> [1] 0.13086102 0.05921165 0.05554814 0.03498161 0.03009382 0.03544119k_tl$K> [1] 0.130861k_tl$P> [1] 2e-04#png("figure4.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure4.pdf")
hist(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "")
abline(v = k_tl$K, lwd = 2, lty = "dotted")
text(x = k_tl$K - 0.02, y = 2500, "Observed value \n of K")
box()
#dev.off()
## PGLS models within classes
## Mammals
pruned_tree_mammals <- drop.tip(pruned_data_tree, pruned_dat$Scientific_name[!pruned_dat$Class == "Mammalia"])
name.check(pruned_tree_mammals, data = pruned_dat[pruned_dat$Class == "Mammalia", ])> [1] "OK"model8 <- gls(new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan, correlation = corBrownian(phy = pruned_tree_mammals, form = ~Scientific_name), data = pruned_dat[pruned_dat$Class == "Mammalia", ], method = "ML")
summary(model8)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan
> Data: pruned_dat[pruned_dat$Class == "Mammalia", ]
> AIC BIC logLik
> 96.58006 107.2957 -43.29003
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.316300 0.8970282 3.696985 0.0005
> log_mass 0.001644 0.0777088 0.021161 0.9832
> log.lifespan 0.003464 0.2564737 0.013504 0.9893
> log_mass:log.lifespan -0.010498 0.0215340 -0.487518 0.6277
>
> Correlation:
> (Intr) lg_mss lg.lfs
> log_mass -0.712
> log.lifespan -0.778 0.682
> log_mass:log.lifespan 0.745 -0.920 -0.853
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.45577217 -0.57238601 -0.26650191 0.08909956 1.42213153
>
> Residual standard error: 0.8803278
> Degrees of freedom: 63 total; 59 residualmodel9 <- update(model8, ~ .-log_mass:log.lifespan)
summary(model9)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan
> Data: pruned_dat[pruned_dat$Class == "Mammalia", ]
> AIC BIC logLik
> 94.83334 103.4059 -43.41667
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.642223 0.5942858 6.128739 0.0000
> log_mass -0.033205 0.0302842 -1.096433 0.2773
> log.lifespan -0.103183 0.1330316 -0.775630 0.4410
>
> Correlation:
> (Intr) lg_mss
> log_mass -0.102
> log.lifespan -0.410 -0.500
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.42754826 -0.56875973 -0.29030500 0.09795673 1.42956182
>
> Residual standard error: 0.8820992
> Degrees of freedom: 63 total; 60 residualmodel10 <- update(model9, ~ .-log.lifespan)
summary(model10)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass
> Data: pruned_dat[pruned_dat$Class == "Mammalia", ]
> AIC BIC logLik
> 93.46188 99.89128 -43.73094
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.453449 0.5403901 6.39066 0.0000
> log_mass -0.044940 0.0261481 -1.71867 0.0907
>
> Correlation:
> (Intr)
> log_mass -0.388
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.4666567 -0.5984019 -0.3109612 0.0536049 1.4274522
>
> Residual standard error: 0.8865104
> Degrees of freedom: 63 total; 61 residual## Birds + Reptiles
bird.reptiles <- rbind(pruned_dat[pruned_dat$Class == "Aves", ], pruned_dat[pruned_dat$Class == "Reptilia", ])
pruned_dat_minus_bird.reptiles <- pruned_dat[!pruned_dat$Scientific_name %in% bird.reptiles$Scientific_name, ]
pruned_tree_bird.reptiles <- drop.tip(pruned_data_tree, pruned_dat_minus_bird.reptiles$Scientific_name)
name.check(pruned_tree_bird.reptiles, bird.reptiles)> [1] "OK"model11 <- gls(new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm, correlation = corBrownian(phy = pruned_tree_bird.reptiles, form = ~Scientific_name), data = bird.reptiles, method = "ML")
summary(model11)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm
> Data: bird.reptiles
> AIC BIC logLik
> 162.1094 179.842 -72.05468
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 4.881621 7.254073 0.6729490 0.5044
> log_mass -0.118243 1.236757 -0.0956076 0.9243
> log.lifespan -0.673792 2.146858 -0.3138505 0.7551
> Endo_ectothermendo -3.252267 8.043313 -0.4043442 0.6879
> log_mass:log.lifespan 0.059966 0.347826 0.1724018 0.8639
> log_mass:Endo_ectothermendo 0.410976 1.317415 0.3119564 0.7565
> log.lifespan:Endo_ectothermendo 1.217913 2.343425 0.5197150 0.6058
> log_mass:log.lifespan:Endo_ectothermendo -0.190072 0.375291 -0.5064671 0.6150
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass -0.900
> log.lifespan -0.951 0.909
> Endo_ectothermendo -0.897 0.807 0.854
> log_mass:log.lifespan 0.871 -0.979 -0.931 -0.778
> log_mass:Endo_ectothermendo 0.845 -0.939 -0.853 -0.873 0.920
> log.lifespan:Endo_ectothermendo 0.871 -0.833 -0.916 -0.922 0.853
> log_mass:log.lifespan:Endo_ectothermendo -0.807 0.908 0.863 0.847 -0.927
> lg_:E_ lg.:E_
> log_mass
> log.lifespan
> Endo_ectothermendo
> log_mass:log.lifespan
> log_mass:Endo_ectothermendo
> log.lifespan:Endo_ectothermendo 0.903
> log_mass:log.lifespan:Endo_ectothermendo -0.977 -0.932
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.01611101 -0.18310917 -0.04275159 0.11800721 0.55779323
>
> Residual standard error: 2.236097
> Degrees of freedom: 53 total; 45 residualmodel12 <- update(model11, ~ .-log_mass:log.lifespan:Endo_ectotherm)
summary(model12)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm
> Data: bird.reptiles
> AIC BIC logLik
> 160.4106 176.1729 -72.2053
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 1.9158767 4.246819 0.4511321 0.6540
> log_mass 0.4503777 0.514506 0.8753592 0.3859
> log.lifespan 0.2646984 1.075332 0.2461552 0.8067
> Endo_ectothermendo 0.1969454 4.244780 0.0463971 0.9632
> log_mass:log.lifespan -0.1033042 0.129554 -0.7973846 0.4293
> log_mass:Endo_ectothermendo -0.2409376 0.278353 -0.8655822 0.3912
> log.lifespan:Endo_ectothermendo 0.1118884 0.843186 0.1326972 0.8950
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:. lg_:E_
> log_mass -0.676
> log.lifespan -0.852 0.593
> Endo_ectothermendo -0.679 0.171 0.459
> log_mass:log.lifespan 0.554 -0.877 -0.693 0.034
> log_mass:Endo_ectothermendo 0.449 -0.580 -0.094 -0.406 0.175
> log.lifespan:Endo_ectothermendo 0.554 0.086 -0.610 -0.691 -0.077 -0.103
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.99688508 -0.16468067 0.01665759 0.12132122 0.57166703
>
> Residual standard error: 2.242461
> Degrees of freedom: 53 total; 46 residualmodel13 <- update(model12, ~ .-log.lifespan:Endo_ectotherm)
summary(model13)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm
> Data: bird.reptiles
> AIC BIC logLik
> 158.4309 172.2229 -72.21545
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 1.6035135 3.497602 0.4584608 0.6487
> log_mass 0.4444743 0.507194 0.8763399 0.3853
> log.lifespan 0.3517726 0.842963 0.4173048 0.6784
> Endo_ectothermendo 0.5859354 3.037761 0.1928840 0.8479
> log_mass:log.lifespan -0.1019766 0.127810 -0.7978772 0.4290
> log_mass:Endo_ectothermendo -0.2371302 0.273962 -0.8655599 0.3911
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass -0.874
> log.lifespan -0.779 0.818
> Endo_ectothermendo -0.492 0.320 0.065
> log_mass:log.lifespan 0.719 -0.876 -0.936 -0.026
> log_mass:Endo_ectothermendo 0.611 -0.576 -0.199 -0.664 0.168
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.00257514 -0.17211069 -0.01640975 0.10795757 0.56656400
>
> Residual standard error: 2.24289
> Degrees of freedom: 53 total; 47 residualmodel14 <- update(model13, ~ .-log_mass:log.lifespan)
summary(model14)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:Endo_ectotherm
> Data: bird.reptiles
> AIC BIC logLik
> 157.144 168.9657 -72.57198
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.610949 2.4204177 1.4918700 0.1423
> log_mass 0.089910 0.2435472 0.3691674 0.7136
> log.lifespan -0.278043 0.2946686 -0.9435777 0.3501
> Endo_ectothermendo 0.522535 3.0252042 0.1727271 0.8636
> log_mass:Endo_ectothermendo -0.200409 0.2690439 -0.7448938 0.4600
>
> Correlation:
> (Intr) lg_mss lg.lfs End_ct
> log_mass -0.727
> log.lifespan -0.433 -0.016
> Endo_ectothermendo -0.682 0.617 0.116
> log_mass:Endo_ectothermendo 0.716 -0.903 -0.120 -0.669
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.98809648 -0.17093673 -0.03260668 0.10129896 0.56675072
>
> Residual standard error: 2.258029
> Degrees of freedom: 53 total; 48 residualmodel15 <- update(model14, ~ .-log_mass:Endo_ectotherm)
summary(model15)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm
> Data: bird.reptiles
> AIC BIC logLik
> 155.7531 165.6046 -72.87655
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 4.902176 1.6815740 2.9152307 0.0053
> log_mass -0.073915 0.1041454 -0.7097328 0.4812
> log.lifespan -0.304342 0.2912141 -1.0450783 0.3011
> Endo_ectothermendo -0.985788 2.2373682 -0.4406018 0.6614
>
> Correlation:
> (Intr) lg_mss lg.lfs
> log_mass -0.267
> log.lifespan -0.500 -0.292
> Endo_ectothermendo -0.390 0.038 0.049
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -0.96251802 -0.14477137 0.01596402 0.15970774 0.58580223
>
> Residual standard error: 2.271042
> Degrees of freedom: 53 total; 49 residualmodel16 <- update(model15, ~ .-Endo_ectotherm)
summary(model16)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan
> Data: bird.reptiles
> AIC BIC logLik
> 153.9627 161.8438 -72.98133
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 4.613051 1.5357249 3.0038264 0.0042
> log_mass -0.072153 0.1032265 -0.6989782 0.4878
> log.lifespan -0.298088 0.2885145 -1.0331819 0.3065
>
> Correlation:
> (Intr) lg_mss
> log_mass -0.274
> log.lifespan -0.523 -0.294
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.28077125 -0.44823059 -0.29285806 -0.05635584 0.50043747
>
> Residual standard error: 2.275537
> Degrees of freedom: 53 total; 50 residualmodel17 <- update(model16, ~ .-log_mass)
summary(model17)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log.lifespan
> Data: bird.reptiles
> AIC BIC logLik
> 152.478 158.3889 -73.23902
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 4.318816 1.4694824 2.939005 0.0049
> log.lifespan -0.357471 0.2743365 -1.303038 0.1984
>
> Correlation:
> (Intr)
> log.lifespan -0.657
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.28381370 -0.36979441 -0.23548939 -0.04838853 0.54703250
>
> Residual standard error: 2.286627
> Degrees of freedom: 53 total; 51 residual## Fish
pruned_tree_fish <- drop.tip(pruned_data_tree, pruned_dat$Scientific_name[!pruned_dat$Class == "Fish"])
name.check(pruned_tree_fish, data = pruned_dat[pruned_dat$Class == "Fish", ])> [1] "OK"model18 <- gls(new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan, correlation = corBrownian(phy = pruned_tree_fish, form = ~Scientific_name), data = pruned_dat[pruned_dat$Class == "Fish", ], method = "ML")
summary(model18)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan
> Data: pruned_dat[pruned_dat$Class == "Fish", ]
> AIC BIC logLik
> 52.46258 58.75307 -21.23129
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 2.8425234 0.6822877 4.166165 0.0004
> log_mass -0.0319071 0.1149743 -0.277515 0.7840
> log.lifespan 0.0608574 0.2242946 0.271328 0.7887
> log_mass:log.lifespan -0.0097255 0.0392377 -0.247860 0.8065
>
> Correlation:
> (Intr) lg_mss lg.lfs
> log_mass -0.420
> log.lifespan -0.568 0.445
> log_mass:log.lifespan 0.421 -0.913 -0.694
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.52656470 -1.07403325 -0.87097331 0.01258059 0.90215128
>
> Residual standard error: 0.9651793
> Degrees of freedom: 26 total; 22 residualmodel19 <- update(model18, ~ .-log_mass:log.lifespan)
summary(model19)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass + log.lifespan
> Data: pruned_dat[pruned_dat$Class == "Fish", ]
> AIC BIC logLik
> 50.53509 55.56747 -21.26754
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 2.9137354 0.6060885 4.807442 0.0001
> log_mass -0.0579381 0.0458255 -1.264320 0.2188
> log.lifespan 0.0222689 0.1581329 0.140824 0.8892
>
> Correlation:
> (Intr) lg_mss
> log_mass -0.094
> log.lifespan -0.423 -0.644
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.5039949921 -1.0703688400 -0.8607380401 -0.0004501222 0.8154978887
>
> Residual standard error: 0.966526
> Degrees of freedom: 26 total; 23 residualmodel20 <- update(model19, ~ .-log.lifespan)
summary(model20)> Generalized least squares fit by maximum likelihood
> Model: new.log.TL ~ log_mass
> Data: pruned_dat[pruned_dat$Class == "Fish", ]
> AIC BIC logLik
> 48.5575 52.33179 -21.27875
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 2.949802 0.5379849 5.483057 0.0000
> log_mass -0.053780 0.0343219 -1.566929 0.1302
>
> Correlation:
> (Intr)
> log_mass -0.529
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.50239115 -1.05995948 -0.87321945 0.01304271 0.83486397
>
> Residual standard error: 0.9669426
> Degrees of freedom: 26 total; 24 residualpruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 0] <- "absent"
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 1] <- "present"
pruned_dat$Telomerase_activity[is.na(pruned_dat$Telomerase_activity)] <- "N/A"
tel.act <- setNames(pruned_dat$Telomerase_activity, rownames(pruned_dat))
tel.act> Pygoscelis_adeliae Tachymarptis_melba
> "N/A" "N/A"
> Haliaeetus_leucocephalus Phoenicopterus_ruber
> "N/A" "N/A"
> Falco_sparverius Ichthyaetus_audouinii
> "N/A" "N/A"
> Hirundo_rustica Branta_leucopsis
> "N/A" "N/A"
> Cepphus_grylle Rissa_tridactyla
> "N/A" "N/A"
> Turdus_merula Gallus_gallus
> "N/A" "absent"
> Sterna_hirundo Calidris_alpina
> "present" "N/A"
> Haematopus_ostralegus Aphelocoma_coerulescens
> "N/A" "N/A"
> Chrysolophus_pictus Fregata_minor
> "N/A" "N/A"
> Parus_major Gyps_fulvus
> "N/A" "N/A"
> Strigops_habroptila Oceanodroma_leucorhoa
> "N/A" "present"
> Larus_fuscus Tachycineta_albilinea
> "N/A" "N/A"
> Colinus_virginianus Fulmarus_glacialis
> "N/A" "N/A"
> Macronectes_halli Accipiter_gentilis
> "N/A" "N/A"
> Amazona_amazonica Buteo_jamaicensis
> "N/A" "N/A"
> Vireo_olivaceus Sula_sula
> "N/A" "N/A"
> Riparia_riparia Passerculus_sandwichensis
> "N/A" "N/A"
> Macronectes_giganteus Buteo_swainsoni
> "N/A" "N/A"
> Uria_lomvia Aphelocoma_ultramarina
> "N/A" "N/A"
> Tachycineta_bicolor Diomedea_exulans
> "present" "N/A"
> Lonchura_striata Meleagris_gallopavo
> "N/A" "N/A"
> Taeniopygia_guttata Calonectris_diomedea
> "present" "N/A"
> Acrocephalus_sechellensis Anguilla_rostrata
> "N/A" "present"
> Oryzias_latipes Pseudocaranx_wrighti
> "present" "N/A"
> Oreochromis_niloticus Carassius_auratus
> "N/A" "N/A"
> Cyprinus_carpio Danio_rerio
> "N/A" "present"
> Xiphophorus_maculatus Xiphophorus_hellerii
> "N/A" "N/A"
> Nothobranchius_furzeri Nothobranchius_rachovii
> "present" "N/A"
> Fundulus_heteroclitus Gadus_morhua
> "present" "present"
> Dicentrarchus_labrax Acanthopagrus_schlegelii
> "N/A" "N/A"
> Upeneichthys_vlamingii Macquaria_ambigua
> "N/A" "N/A"
> Lutjanus_argentimaculatus Oncorhynchus_mykiss
> "N/A" "present"
> Platycephalus_bassensis Mustelus_asterias
> "N/A" "N/A"
> Callorhinchus_milii Urolophus_paucimaculatus
> "N/A" "N/A"
> Squalus_acanthias Torpedo_marmorata
> "present" "N/A"
> Neoceratodus_forsteri Setifer_setosus
> "N/A" "present"
> Balaena_mysticetus Eschrichtius_robustus
> "absent" "absent"
> Giraffa_camelopardalis Ovis_aries
> "absent" "absent"
> Rangifer_tarandus Capra_hircus
> "N/A" "N/A"
> Tursiops_truncatus Camelus_dromedarius
> "absent" "present"
> Sus_scrofa Muntiacus_reevesi
> "present" "absent"
> Muntiacus_muntjak Bos_taurus
> "absent" "absent"
> Zalophus_californianus Crocuta_crocuta
> "absent" "absent"
> Canis_lupus Ursus_maritimus
> "absent" "absent"
> Felis_catus Ailurus_fulgens
> "absent" "absent"
> Panthera_tigris Meles_meles
> "present" "N/A"
> Myotis_lucifugus Tadarida_brasiliensis
> "present" "present"
> Pteropus_rodricensis Chaetophractus_vellerosus
> "absent" "absent"
> Didelphis_virginiana Atelerix_albiventris
> "present" "present"
> Sorex_araneus Procavia_capensis
> "present" "absent"
> Lepus_californicus Sylvilagus_aquaticus
> "absent" "absent"
> Oryctolagus_cuniculus Ochotona_princeps
> "absent" "present"
> Macroscelides_proboscideus Galegeeska_rufescens
> "present" "present"
> Ceratotherium_simum Equus_grevyi
> "absent" "absent"
> Equus_caballus Tapirus_indicus
> "absent" "absent"
> Myrmecophaga_tridactyla Choloepus_hoffmanni
> "absent" "absent"
> Homo_sapiens Ateles_geoffroyi
> "absent" "absent"
> Saimiri_sciureus Pan_paniscus
> "absent" "absent"
> Pongo_pygmaeus Pan_troglodytes
> "absent" "N/A"
> Gorilla_gorilla Lemur_catta
> "N/A" "absent"
> Macaca_nemestrina Macaca_mulatta
> "N/A" "absent"
> Macaca_fascicularis Loxodonta_africana
> "absent" "absent"
> Elephas_maximus Castor_canadensis
> "absent" "absent"
> Hydrochoerus_hydrochaeris Sciurus_carolinensis
> "absent" "present"
> Aplodontia_rufa Heterocephalus_glaber
> "absent" "present"
> Rattus_norvegicus Mus_musculus
> "present" "present"
> Tupaia_tana Tupaia_belangeri
> "present" "N/A"
> Alligator_sinensis Alligator_mississippiensis
> "absent" "N/A"
> Zootoca_vivipara Lacerta_agilis
> "N/A" "N/A"
> Liasis_fuscus Emys_orbicularis
> "N/A" "N/A"
> Trachemys_scripta Chrysemys_picta
> "N/A" "present"TA <- to.matrix(tel.act, unique(tel.act))
TA <- TA[pruned_data_tree$tip.label, ]
life.span <- setNames(pruned_dat$log.lifespan, rownames(pruned_dat))
log_mass <- setNames(pruned_dat$log_mass, rownames(pruned_dat))
telo.length <- setNames(pruned_dat$new.log.TL, rownames(pruned_dat))plotTree(pruned_data_tree, ftype = "off", mar = c(3, 2, 2, 3))
tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.22, offset = 4.3)
par(xpd = TRUE)
legend(x = 150, y = 0.2, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)
par(new = TRUE)
par(mar = c(3, 32, 2, 1.1))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1, space = 0,
ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)
axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)
Figure 3
## Ancestral state reconstruction of telomere size
fit <- fastAnc(pruned_data_tree, telo.length, vars = TRUE, CI = TRUE)
fit$CI[1, ]> [1] 0.9673174 4.3349363obj <- contMap(pruned_data_tree, telo.length, plot = FALSE)
##png("figure5.png", width = 7, height = 7, units = "in", res = 360)
##pdf("figure5.pdf")
plot(obj, ftype = "off", legend = FALSE, ylim = c(1-0.09*(Ntip(obj$tree)-1), Ntip(obj$tree)), mar = c(1, 0.1, 1, 5), lwd = 1.5)
add.color.bar(150, obj$cols,title = "Log telomere length (kb)", lims = obj$lims, digits = 3, prompt = FALSE, x = 0,
y = 1-0.08*(Ntip(obj$tree)-1), lwd = 4, fsize = 0.6, subtitle = "")
tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.16, offset = 5.7)
legend(x = 200, y = -5, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)
par(new = TRUE)
par(mar = c(3.6, 29.7, 3, 3))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)
axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)
par(new = TRUE)
par(mar = c(3.6, 32.2, 3, 0.5))
barplot(log_mass[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)
axis(1, at = round(seq(min(log_mass), max(log_mass), 5), 1), labels = FALSE)
text(round(seq(min(log_mass), max(log_mass), 5), 1), par("usr")[3] - 0.2, labels = round(seq(min(log_mass), max(log_mass), 5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log mass \n (gr)", side = 1, line = 1.6, cex = 0.5, font = 2)
##dev.off()## Correlated evolution under the threshold model
## Removing NAs from the dataset
pruned_dat$Telomerase_activity> [1] "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A"
> [8] "N/A" "N/A" "N/A" "N/A" "absent" "present" "N/A"
> [15] "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A"
> [22] "present" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A"
> [29] "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A"
> [36] "N/A" "N/A" "N/A" "present" "N/A" "N/A" "N/A"
> [43] "present" "N/A" "N/A" "present" "present" "N/A" "N/A"
> [50] "N/A" "N/A" "present" "N/A" "N/A" "present" "N/A"
> [57] "present" "present" "N/A" "N/A" "N/A" "N/A" "N/A"
> [64] "present" "N/A" "N/A" "N/A" "N/A" "present" "N/A"
> [71] "N/A" "present" "absent" "absent" "absent" "absent" "N/A"
> [78] "N/A" "absent" "present" "present" "absent" "absent" "absent"
> [85] "absent" "absent" "absent" "absent" "absent" "absent" "present"
> [92] "N/A" "present" "present" "absent" "absent" "present" "present"
> [99] "present" "absent" "absent" "absent" "absent" "present" "present"
> [106] "present" "absent" "absent" "absent" "absent" "absent" "absent"
> [113] "absent" "absent" "absent" "absent" "absent" "N/A" "N/A"
> [120] "absent" "N/A" "absent" "absent" "absent" "absent" "absent"
> [127] "absent" "present" "absent" "present" "present" "present" "present"
> [134] "N/A" "absent" "N/A" "N/A" "N/A" "N/A" "N/A"
> [141] "N/A" "present"pruned_dat_not_NAs <- pruned_dat[!pruned_dat$Telomerase_activity == "N/A", ]
pruned_dat_not_NAs$Telomerase_activity> [1] "absent" "present" "present" "present" "present" "present" "present"
> [8] "present" "present" "present" "present" "present" "present" "present"
> [15] "absent" "absent" "absent" "absent" "absent" "present" "present"
> [22] "absent" "absent" "absent" "absent" "absent" "absent" "absent"
> [29] "absent" "absent" "present" "present" "present" "absent" "absent"
> [36] "present" "present" "present" "absent" "absent" "absent" "absent"
> [43] "present" "present" "present" "absent" "absent" "absent" "absent"
> [50] "absent" "absent" "absent" "absent" "absent" "absent" "absent"
> [57] "absent" "absent" "absent" "absent" "absent" "absent" "absent"
> [64] "present" "absent" "present" "present" "present" "present" "absent"
> [71] "present"check2 <- name.check(pruned_data_tree, pruned_dat_not_NAs)
rm_phy2 <- check2$tree_not_data
rm_dat2 <- check2$data_not_tree
pruned_data_tree_not_NAs <- drop.tip(pruned_data_tree, rm_phy2)
pruned_dat_not_NAs <- subset(pruned_dat_not_NAs, subset = pruned_dat_not_NAs$Scientific_name %in% pruned_data_tree_not_NAs$tip, select = names(pruned_dat_not_NAs))
pruned_dat_not_NAs$Telomerase_activity <- as.factor(pruned_dat_not_NAs$Telomerase_activity)
str(pruned_dat_not_NAs)> 'data.frame': 71 obs. of 20 variables:
> $ Common.name : chr "Chicken " "common tern " "Leach storm petrel " "tree swallow " ...
> $ Scientific_name : chr "Gallus_gallus" "Sterna_hirundo" "Oceanodroma_leucorhoa" "Tachycineta_bicolor" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Aves" "Aves" "Aves" "Aves" ...
> $ Order : chr "Galliformes" "Charadriiformes" "Procellariiformes" "Passeriformes" ...
> $ Family : chr "Phasianidae" "Sternidae" "Hydrobatidae" "Hirundinidae" ...
> $ Genus : chr "Gallus" "Sterna" "Oceanodroma" "Tachynicineta" ...
> $ Species : chr "gallus" "hirundo" "leucorhoa" "bicolor" ...
> $ Endo_ectotherm : chr "endo" "endo" "endo" "endo" ...
> $ Adult_mass_grams : num 779.8 120 44.6 19 12 ...
> $ log_mass : num 6.66 4.8 3.82 3 2.56 ...
> $ Lifespan_years : num 30 33 36 12.1 12 50 5 5.5 0.23 4 ...
> $ Loglifespan : num 1.48 1.52 1.56 1.08 1.08 ...
> $ Average_Telomere_Length_kb: num 19.7 12.8 19 13.6 13.7 ...
> $ logTL : num 1.29 1.11 1.28 1.13 1.14 ...
> $ Telomerase_activity : Factor w/ 2 levels "absent","present": 1 2 2 2 2 2 2 2 2 2 ...
> $ log.lifespan : num 3.43 3.53 3.61 2.57 2.56 ...
> $ new.log.TL : num 2.98 2.55 2.95 2.61 2.62 ...head(pruned_dat_not_NAs)> Common.name Scientific_name Domain
> Gallus_gallus Chicken Gallus_gallus Eukaryota
> Sterna_hirundo common tern Sterna_hirundo Eukaryota
> Oceanodroma_leucorhoa Leach storm petrel Oceanodroma_leucorhoa Eukaryota
> Tachycineta_bicolor tree swallow Tachycineta_bicolor Eukaryota
> Taeniopygia_guttata Zebra finch Taeniopygia_guttata Eukaryota
> Anguilla_rostrata American eel Anguilla_rostrata Eukaryota
> Kingdom Phylum Class Order Family
> Gallus_gallus Metazoa Chordata Aves Galliformes Phasianidae
> Sterna_hirundo Metazoa Chordata Aves Charadriiformes Sternidae
> Oceanodroma_leucorhoa Metazoa Chordata Aves Procellariiformes Hydrobatidae
> Tachycineta_bicolor Metazoa Chordata Aves Passeriformes Hirundinidae
> Taeniopygia_guttata Metazoa Chordata Aves Passeriformes Estrildidae
> Anguilla_rostrata Metazoa Chordata Fish Anguilliformes Anguilla
> Genus Species Endo_ectotherm Adult_mass_grams
> Gallus_gallus Gallus gallus endo 779.8
> Sterna_hirundo Sterna hirundo endo 120.0
> Oceanodroma_leucorhoa Oceanodroma leucorhoa endo 44.6
> Tachycineta_bicolor Tachynicineta bicolor endo 19.0
> Taeniopygia_guttata Taeniopygia guttata endo 12.0
> Anguilla_rostrata Anguilla rostrata ecto 4032.0
> log_mass Lifespan_years Loglifespan
> Gallus_gallus 6.660319 30.0 1.477121
> Sterna_hirundo 4.795791 33.0 1.518514
> Oceanodroma_leucorhoa 3.819908 36.0 1.556303
> Tachycineta_bicolor 2.995732 12.1 1.082785
> Taeniopygia_guttata 2.564949 12.0 1.079181
> Anguilla_rostrata 8.302266 50.0 1.698970
> Average_Telomere_Length_kb logTL Telomerase_activity
> Gallus_gallus 19.70 1.294466 absent
> Sterna_hirundo 12.80 1.107210 present
> Oceanodroma_leucorhoa 19.03 1.279439 present
> Tachycineta_bicolor 13.60 1.133539 present
> Taeniopygia_guttata 13.70 1.136721 present
> Anguilla_rostrata 12.50 1.096910 present
> log.lifespan new.log.TL
> Gallus_gallus 3.433987 2.980619
> Sterna_hirundo 3.526361 2.549445
> Oceanodroma_leucorhoa 3.610918 2.946017
> Tachycineta_bicolor 2.572612 2.610070
> Taeniopygia_guttata 2.564949 2.617396
> Anguilla_rostrata 3.931826 2.525729pruned_data_tree_not_NAs>
> Phylogenetic tree with 71 tips and 70 internal nodes.
>
> Tip labels:
> Squalus_acanthias, Oncorhynchus_mykiss, Oryzias_latipes, Fundulus_heteroclitus, Nothobranchius_furzeri, Gadus_morhua, ...
> Node labels:
> , '286', '79', '80', '87', '84', ...
>
> Rooted; includes branch lengths.name.check(pruned_data_tree_not_NAs, pruned_dat_not_NAs)> [1] "OK"names(pruned_dat_not_NAs)> [1] "Common.name" "Scientific_name"
> [3] "Domain" "Kingdom"
> [5] "Phylum" "Class"
> [7] "Order" "Family"
> [9] "Genus" "Species"
> [11] "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years"
> [15] "Loglifespan" "Average_Telomere_Length_kb"
> [17] "logTL" "Telomerase_activity"
> [19] "log.lifespan" "new.log.TL"## Set the number of generations
ngen <- 5e6
## Run MCMC
mcmc.model <- threshBayes(pruned_data_tree_not_NAs, pruned_dat_not_NAs[ , c(17, 18)], type = c("cont", "disc"), ngen = ngen,
plot = FALSE, control = list(print.interval = 5e+05))> Starting MCMC....
> generation: 500000; mean acceptance rate: 0.26
> generation: 1000000; mean acceptance rate: 0.24
> generation: 1500000; mean acceptance rate: 0.26
> generation: 2000000; mean acceptance rate: 0.22
> generation: 2500000; mean acceptance rate: 0.26
> generation: 3000000; mean acceptance rate: 0.24
> generation: 3500000; mean acceptance rate: 0.22
> generation: 4000000; mean acceptance rate: 0.24
> generation: 4500000; mean acceptance rate: 0.25
> generation: 5000000; mean acceptance rate: 0.22
> Done MCMC.mcmc.model>
> Object of class "threshBayes" consisting of a matrix (L) of
> sampled liabilities for the tips of the tree & a second matrix
> (par) with the sample model parameters & correlation.
>
> Mean correlation (r) from the posterior sample is: 0.478.
>
> Ordination of discrete traits:
>
> Trait 2: absent <-> present## Pull out the post burn-in sample and compute HPD
r.mcmc <- tail(mcmc.model$par$r, 0.8*nrow(mcmc.model$par))
class(r.mcmc) <- "mcmc"
hpd.r <- HPDinterval(r.mcmc)
hpd.r> lower upper
> var1 0.1810489 0.7420291
> attr(,"Probability")
> [1] 0.9500012
Figure S1
## Profile plots from a Bayesian MCMC analysis of the threshold model
#png("figureS1.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS1.pdf")
plot(mcmc.model)
#dev.off()
Figure 6
## Plot posterior density
#png("figure6.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure6.pdf")
par(las = 1)
plot(density(mcmc.model), xlim = c(-1, 1.5))
## add whiskers to show HPD
h <- 0-par()$usr[3]
lines(x = hpd.r, y = rep(-h/2, 2))
lines(x = rep(hpd.r[1], 2), y = c(-0.3, -0.7)*h)
lines(x = rep(hpd.r[2], 2), y = c(-0.3, -0.7)*h)
box()
##dev.off()